Take_home_Ex1

Overview

The recent shift in payment being made more digital, companies and organisations can more easily collect data and information that are linked to consumer habits. The transportation industry including public transport such as buses has also lean into this phenomenon. The information collected include travelling patterns that can help companies plan for more efficient routes or where heavy ridership is to be expected.

Objectives

Exploratory Spatial Data Analysis (ESDA) hold tremendous potential to address complex problems facing society. In this study, you are tasked to apply appropriate Local Indicators of Spatial Association (GLISA) and Emerging Hot Spot Analysis (EHSA) to undercover the spatial and spatio-temporal mobility patterns of public bus passengers in Singapore.

Task

Here we will utilise bus travelling data at different time duration for plotting out geospatial data and analysing them using various statistical tools.

Geovisualisation and Analysis

Computing passenger trips at the hexagonal level for the following time intervals:

• Weekday morning peak, 6am to 9am
• Weekday afternoon peak, 5pm to 8pm
• Weekend/holiday morning peak, 11am to 2pm
• Weekend/holiday evening peak, 4pm to 7pm

Display the geographical distribution using choropleth maps of the hexagons.

Combine all of the passenger trips made by all of the bus stops within a hexagon together

Local Indicators of Spatial Association(LISA) Analysis

Utilise Queen’s contiguity for performing LISA of the passenger trips by origin at hexagonal level Displat the LISA maps of the passenger trips at hexagonal level.

Load Packages and Data

Load packages

Here we will load the packages needed for this exercise and their respective functions - sf: - tmap: - spdep: - tidyverse: - dplyr: - mapview: - sfdep:

pacman::p_load(sf,tmap,spdep,tidyverse, dplyr, mapview, sfdep)

Loading data

Loading aspatial table

Here we will read all of the ridership from different bus stops in Oct 2023 and assign it to the variable.

odbus <- read_csv("data/aspatial/origin_destination_bus_202310.csv")

We will then extract the information from the following and assign them to different variables.

Day	Duration	Variable name
Weekdays	6am - 9am	weekdayAM_6_9
Weekdays	5pm - 8pm	weekdayPM_5_8
Weekends/Holidays	11am - 2pm	weekendAM_11_14
Weekends/Holidays	4pm - 7pm	weekendPM_4_7

# Filter data for weekday morning hours
weekdayAM_6_9 <- odbus %>% 
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 6 & TIME_PER_HOUR <= 9) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))

weekdayPM_5_8 <- odbus %>% 
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 17 & TIME_PER_HOUR <= 20) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))

weekendAM_11_14 <- odbus %>% 
  filter(DAY_TYPE == "WEEKENDS/HOLIDAY") %>%
  filter(TIME_PER_HOUR >= 11 & TIME_PER_HOUR <= 14) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))

weekendPM_16_19 <- odbus %>% 
  filter(DAY_TYPE == "WEEKENDS/HOLIDAY") %>%
  filter(TIME_PER_HOUR >= 17 & TIME_PER_HOUR <= 19) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))

Loading Geospatial data

Next we will import all of the bus stops and their coordinates and attached it to the busstop variable.

# Import geospatial data
busstop <- st_read(dsn = "data/geospatial", layer = "BusStop") %>%
  st_transform(crs = 3414)

Reading layer `BusStop' from data source 
  `C:\Users\Lian Khye\Desktop\MITB\Geospatial\geartooth\ISSS624\Take_home_Ex1\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 5161 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21

We will first rename the bus stop column title for easier data joining.

colnames(busstop)[colnames(busstop) == "BUS_STOP_N"] <- "ORIGIN_PT_CODE"

After that we will create the hexagons that will create the map layout. The hexagons will be shaped 250 x 250 cell size. All of the hexagons will also be given a grid id name that can be used for identifying each individual grid.

center <- st_centroid(busstop)

area_honeycomb_grid <- st_make_grid(center, cellsize = c(250 * sqrt(3), 250 * 2), what = "polygons", square = FALSE)
honeycomb_grid_sf <- st_sf(area_honeycomb_grid) %>%
  mutate(grid_id = 1:length(lengths(area_honeycomb_grid)))

Data processing

Assigning individual bus stop to hexagons

First we will assign the bus stop point geometry data to each polygon using st_intersection(). The function assigns all of the points to a polygon and then join both tables together.

busstop_hex <- st_intersection(busstop, honeycomb_grid_sf) %>%
  st_drop_geometry()

Duplication check

Here we will check for the presence of any duplication before we further process the data.

duplicate1 <- weekdayAM_6_9 %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

duplicate2 <- weekdayPM_5_8 %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

duplicate3 <- weekendAM_11_14 %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

duplicate4 <- weekendPM_16_19 %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

We can see which data points are duplicated here.

c(duplicate1,duplicate2,duplicate3,duplicate4)

$ORIGIN_PT_CODE
character(0)

$TRIPS
numeric(0)

$ORIGIN_PT_CODE
character(0)

$TRIPS
numeric(0)

$ORIGIN_PT_CODE
character(0)

$TRIPS
numeric(0)

$ORIGIN_PT_CODE
character(0)

$TRIPS
numeric(0)

Finally we only keep data points that are unique using the unique() function.

unique_weekdayAM <- unique(weekdayAM_6_9)
unique_weekdayPM <- unique(weekdayPM_5_8)
unique_weekendAM <- unique(weekendAM_11_14)
unique_weekendPM <- unique(weekendPM_16_19)

Join tables

Next we will then join the variables that we created earlier that contains the total number of trips at different time intervals and the busstop_hex variable together using BUS_STOP_N column title that we have in common.

count_weekdayAM_6_9 <- left_join(unique_weekdayAM , busstop_hex) %>%
  group_by(grid_id) %>%
  summarise(TOT_TRIPS = sum(TRIPS))

count_weekdayPM_5_8 <- left_join(unique_weekdayPM , busstop_hex) %>%
  group_by(grid_id) %>%
  summarise(TOT_TRIPS = sum(TRIPS))

count_weekendAM_11_14 <- left_join(unique_weekendAM , busstop_hex) %>%
  group_by(grid_id) %>%
  summarise(TOT_TRIPS = sum(TRIPS))

count_weekendPM_16_19 <- left_join(unique_weekendPM , busstop_hex) %>%
  group_by(grid_id) %>%
  summarise(TOT_TRIPS = sum(TRIPS))

poly_weekdayAM_6_9 <- left_join(honeycomb_grid_sf,count_weekdayAM_6_9)
poly_weekdayPM_5_8 <- left_join(honeycomb_grid_sf,count_weekdayPM_5_8)
poly_weekendAM_11_14 <- left_join(honeycomb_grid_sf,count_weekendAM_11_14)
poly_weekendPM_16_19 <- left_join(honeycomb_grid_sf,count_weekendPM_16_19)

grid_weekdayAM <- poly_weekdayAM_6_9 %>%
  filter(TOT_TRIPS > 0)

grid_weekdayPM <- poly_weekdayPM_5_8 %>%
  filter(TOT_TRIPS > 0)

grid_weekendAM <- poly_weekendAM_11_14 %>%
  filter(TOT_TRIPS > 0)

grid_weekendPM <- poly_weekendPM_16_19 %>%
  filter(TOT_TRIPS > 0)

Choropleth map

Here we will plot the choropleth map for the different time intervals. We will be using tmap_mode(“plot”) to create an interactive map. Although we will be coding in accessories such as the compass, they will not be displayed in the interactive map. However by writing them first, we can display them in subsequent maps once we view them in plot mode.

Weekday 6am to 9am

tmap_mode("view")
mapA <- tm_shape(grid_weekdayAM)+
  tm_fill("TOT_TRIPS", 
          style = "quantile", 
          palette = "Blues",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekday 6am-9am",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))
mapA

The total number of ridership range from 1 to 357043 per hexagon. The range of the data are divided to quantile range bands for clearer distinction between ridership of each hexagon.

Weekday 5pm to 8pm

mapB <- tm_shape(grid_weekdayPM)+
  tm_fill("TOT_TRIPS", 
          style = "quantile", 
          palette = "Reds",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekday 5pm-8pm",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))
mapB

Weekend/Holidays 11am to 2pm

mapC <- tm_shape(grid_weekendAM)+
  tm_fill("TOT_TRIPS", 
          style = "quantile", 
          palette = "Greens",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekend/holidays 11am-2pm",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))
mapC

Weekend/Holidays 4pm to 7pm

mapD <- tm_shape(grid_weekendPM)+
  tm_fill("TOT_TRIPS", 
          style = "quantile", 
          palette = "Purples",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekend/holidays 4pm-7pm",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))
mapD

Weekday 6am to 9am

tmap_mode("view")
mapA <- tm_shape(grid_weekdayAM)+
  tm_fill("TOT_TRIPS", 
          style = "quantile", 
          palette = "Blues",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekday 6am-9am",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))
mapA

The total number of ridership range from 1 to 357043 per hexagon. The range of the data are divided to quantile range bands for clearer distinction between ridership of each hexagon.

Weekday 5pm to 8pm

mapB <- tm_shape(grid_weekdayPM)+
  tm_fill("TOT_TRIPS", 
          style = "quantile", 
          palette = "Reds",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekday 5pm-8pm",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))
mapB

The total number of ridership range from 1 to 568845 per hexagon.

Weekend/Holidays 11am to 2pm

mapC <- tm_shape(grid_weekendAM)+
  tm_fill("TOT_TRIPS", 
          style = "quantile", 
          palette = "Greens",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekend/holidays 11am-2pm",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))
mapC

The total number of ridership range from 1 to 117609 per hexagon.

Weekend/Holidays 4pm to 7pm

mapD <- tm_shape(grid_weekendPM)+
  tm_fill("TOT_TRIPS", 
          style = "quantile", 
          palette = "Purples",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekend/holidays 4pm-7pm",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))
mapD

The total number of ridership range from 1 to 114410 per hexagon.

Plot maps

We can also utilise a plot mode by using tmap_mode(“plot”).

tmap_mode("plot")
mapA

This allows for other accessories to be added such as compass and scales which might be useful depending on the application or use of the map.

Choropleth map discussion

We can see that the total ridership over the weekends is lower than the weekday counterparts despite both being classified as peak hours. This difference is likely due to people travelling to or from work on the weekdays as compared to the weekends where there will be less traffic as people choose to stay at home or electing to travel to districts where it may be more accessible by other transportations such as trains.

LISA

Before we calculate LISA, we will add the number of neighbours to each dataset for the various grids. We will use st_contiguity() and add a neighbour column to the data points under a new variable name. By default, the Queen method of contiguity will be used for calculating the neighbours.

wm_qA <- grid_weekdayAM %>%
    mutate(nb = st_contiguity(area_honeycomb_grid),
         .before = 1)

wm_qB <- grid_weekdayPM %>%
  mutate(nb = st_contiguity(area_honeycomb_grid),
         .before = 1)

wm_qC <- grid_weekendAM %>%
  mutate(nb = st_contiguity(area_honeycomb_grid),
         .before = 1)

wm_qD <- grid_weekendPM %>%
  mutate(nb = st_contiguity(area_honeycomb_grid),
         .before = 1)

Here we can use the summary() function to take a look at the different regions and the neighbours that they have.

summary(wm_qA$nb)

Neighbour list object:
Number of regions: 1779 
Number of nonzero links: 7672 
Percentage nonzero weights: 0.2424134 
Average number of links: 4.312535 
14 regions with no links:
1 311 358 489 526 651 765 837 1226 1753 1756 1760 1769 1779
Link number distribution:

  0   1   2   3   4   5   6 
 14  54 165 250 392 454 450 
54 least connected regions:
2 3 8 22 23 35 36 56 65 70 110 126 172 199 200 209 210 218 235 251 347 373 652 664 698 735 756 806 818 819 844 846 861 898 927 983 1168 1180 1181 1183 1211 1261 1263 1324 1434 1497 1506 1523 1542 1554 1629 1663 1761 1762 with 1 link
450 most connected regions:
29 48 71 78 81 82 87 89 94 95 112 119 124 125 131 135 141 144 148 178 185 186 187 194 195 196 204 205 206 216 217 224 230 241 247 249 256 258 259 265 269 271 272 276 277 278 283 284 285 316 320 331 337 339 340 341 344 349 350 353 354 355 356 362 365 366 369 381 382 383 385 386 397 398 401 408 410 411 427 436 440 457 465 471 474 481 486 487 496 502 503 504 510 514 515 516 522 528 530 531 539 546 558 559 568 572 579 580 583 584 589 590 601 612 613 615 621 624 625 626 631 632 635 636 641 646 647 654 655 659 666 667 670 671 673 674 683 684 687 690 691 700 701 705 710 711 712 717 724 725 730 731 737 742 743 748 752 761 762 766 771 780 781 786 792 797 798 801 809 810 813 821 822 826 827 835 840 849 855 866 867 870 872 882 883 889 900 901 903 907 910 912 916 917 932 936 937 938 952 953 957 958 959 965 970 976 977 978 979 980 981 984 985 991 997 1000 1001 1002 1005 1006 1012 1015 1016 1019 1020 1022 1023 1024 1025 1031 1032 1036 1037 1038 1039 1041 1042 1043 1048 1057 1058 1061 1062 1063 1066 1067 1071 1076 1080 1081 1088 1089 1090 1092 1096 1097 1098 1100 1103 1104 1105 1111 1114 1115 1118 1121 1131 1132 1134 1136 1145 1146 1147 1148 1159 1160 1161 1172 1173 1179 1185 1186 1187 1204 1206 1207 1209 1221 1222 1223 1229 1230 1231 1232 1233 1239 1240 1244 1245 1246 1251 1255 1257 1258 1271 1274 1282 1283 1284 1287 1289 1290 1296 1297 1298 1303 1305 1311 1316 1317 1318 1322 1329 1330 1335 1337 1338 1339 1343 1349 1350 1356 1361 1363 1366 1368 1369 1371 1375 1382 1383 1384 1385 1386 1387 1391 1394 1395 1396 1397 1398 1403 1405 1408 1409 1410 1411 1418 1419 1422 1423 1424 1425 1426 1427 1429 1431 1437 1438 1439 1447 1451 1452 1453 1461 1465 1470 1484 1487 1495 1498 1502 1508 1513 1516 1536 1537 1544 1545 1549 1551 1558 1559 1564 1565 1573 1574 1582 1583 1585 1586 1587 1591 1594 1595 1598 1599 1601 1602 1608 1611 1612 1617 1620 1621 1622 1623 1627 1635 1644 1651 1656 1658 1659 1667 1668 1670 1671 1674 1681 1683 1684 1692 1694 1696 1698 1700 1704 1705 1706 1708 1710 1711 1715 1718 1724 1725 1728 1729 1732 with 6 links

We see that there are 1779 data points and they have an average of 4.3 neighbours. There are 14 hexagons with 0 neighbours and 450 hexagons with 6 neighbours.

Next we will have to remove all of the data points with 0 neighbours using filter and purrr package.

wm_qA_process <- wm_qA %>%
   filter(!purrr::map_lgl(nb, ~ all(. == 0)))
wm_qB_process <- wm_qB %>%
   filter(!purrr::map_lgl(nb, ~ all(. == 0)))
wm_qC_process <- wm_qC %>%
   filter(!purrr::map_lgl(nb, ~ all(. == 0)))
wm_qD_process <- wm_qD %>%
   filter(!purrr::map_lgl(nb, ~ all(. == 0)))

Spatial weights

wm_qA_weighted <- wm_qA_process %>% 
  mutate(nb = st_contiguity(area_honeycomb_grid),
         wt = st_weights(nb, style = "W"),
         .before = 1)
wm_qB_weighted <- wm_qB_process %>% 
  mutate(nb = st_contiguity(area_honeycomb_grid),
         wt = st_weights(nb, style = "W"),
         .before = 1)
wm_qC_weighted <- wm_qC_process %>% 
  mutate(nb = st_contiguity(area_honeycomb_grid),
         wt = st_weights(nb, style = "W"),
         .before = 1)
wm_qD_weighted <- wm_qD_process %>% 
  mutate(nb = st_contiguity(area_honeycomb_grid),
         wt = st_weights(nb, style = "W"),
         .before = 1)

Local Moran’s I

Here we will be using the local_moran() function to calculate the local Moran’s I for each region or county.

lisaA <- wm_qA_weighted %>% 
  mutate(local_moran = local_moran(
    TOT_TRIPS, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_moran)
lisaB <- wm_qB_weighted %>% 
  mutate(local_moran = local_moran(
    TOT_TRIPS, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_moran)
lisaC <- wm_qC_weighted %>% 
  mutate(local_moran = local_moran(
    TOT_TRIPS, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_moran)
lisaD <- wm_qD_weighted %>% 
  mutate(local_moran = local_moran(
    TOT_TRIPS, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_moran)

The output will be a data fram containing the ii, eii, var_ii, z_ii, p_ii, p_ii_sim and p_folded_sum.

Combined visualisation of local Moran’s I and p-value

Here we will place the maps next to each other.

Weekday 6am to 9am

tmap_mode("plot")
map1A <- tm_shape(lisaA) +
  tm_fill("ii") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "local Moran's I of the total number of trips",
            main.title.size = 0.8)

map2A <- tm_shape(lisaA) +
  tm_fill("p_ii",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig")) + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of local Moran's I",
            main.title.size = 0.8)

tmap_arrange(map1A, map2A, ncol = 2)

Weekday 5pm to 8pm

tmap_mode("plot")
map1B <- tm_shape(lisaB) +
  tm_fill("ii") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "local Moran's I of the total number of trips",
            main.title.size = 0.8)

map2B <- tm_shape(lisaB) +
  tm_fill("p_ii",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig")) + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of local Moran's I",
            main.title.size = 0.8)

tmap_arrange(map1B, map2B, ncol = 2)

Weekday 6am to 9am

tmap_mode("plot")
map1C <- tm_shape(lisaC) +
  tm_fill("ii") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "local Moran's I of the total number of trips",
            main.title.size = 0.8)

map2C <- tm_shape(lisaC) +
  tm_fill("p_ii",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig")) + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of local Moran's I",
            main.title.size = 0.8)

tmap_arrange(map1C, map2C, ncol = 2)

Weekday 6am to 9am

tmap_mode("plot")
map1D <- tm_shape(lisaD) +
  tm_fill("ii") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "local Moran's I of the total number of trips",
            main.title.size = 0.8)

map2D <- tm_shape(lisaD) +
  tm_fill("p_ii",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig")) + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of local Moran's I",
            main.title.size = 0.8)

tmap_arrange(map1D, map2D, ncol = 2)

Visualisation of LISA

Here will will visualise LISA where we can see the presence of outliers and clusters. More information can be found here. The following is a newer method for calculating LISA, and require shorter and more concise steps such as not having to manually form the high-high, high-low, low-high and low-low quadrants. Just make sure that if the data is skewed, we will have to use the median for forming the quadrant.

Weekday 6am to 9am

lisa_sigA <- lisaA  %>%
  filter(p_ii < 0.05)
tmap_mode("view")
tm_shape(lisaA) +
  tm_polygons() +
  tm_borders(alpha = 0.5) +
tm_shape(lisa_sigA) +
  tm_fill("mean") + 
  tm_borders(alpha = 0.4)

Weekday 5pm to 9pm

lisa_sigB <- lisaB  %>%
  filter(p_ii < 0.05)
tmap_mode("view")
tm_shape(lisaB) +
  tm_polygons() +
  tm_borders(alpha = 0.5) +
tm_shape(lisa_sigB) +
  tm_fill("mean") + 
  tm_borders(alpha = 0.4)

Weekends/Holidays 11am to 2pm

lisa_sigC <- lisaC  %>%
  filter(p_ii < 0.05)
tmap_mode("view")
tm_shape(lisaC) +
  tm_polygons() +
  tm_borders(alpha = 0.5) +
tm_shape(lisa_sigC) +
  tm_fill("mean") + 
  tm_borders(alpha = 0.4)

Weekends/Holidays 11am to 2pm

lisa_sigD <- lisaD  %>%
  filter(p_ii < 0.05)
tmap_mode("view")
tm_shape(lisaD) +
  tm_polygons() +
  tm_borders(alpha = 0.5) +
tm_shape(lisa_sigD) +
  tm_fill("mean") + 
  tm_borders(alpha = 0.4)

Checking consistency

lisaA_test <- lisaA %>%
  select(10, 15) %>%
  st_drop_geometry()
colnames(lisaA_test)[colnames(lisaA_test) == "mean"] <- "Weekday AM"


lisaC_test <- lisaC %>%
  select(10, 15) %>%
  st_drop_geometry()
colnames(lisaC_test)[colnames(lisaC_test) == "mean"] <- "Weekend AM"

AM <-  left_join(lisaA_test , lisaC_test) 
AM <- mutate(AM, is_equal = `Weekday AM` == `Weekend AM`)

colnames(AM)[4] = "equality"

AM_diff<-AM$grid_id[AM[4] == FALSE]
AM_diff_unique<-unique(AM_diff)
p <- data.frame(AM_diff_unique)

#origin_SZ <- left_join(origin7_9 , busstop_mpsz,
#            by = c("ORIGIN_PT_CODE" = "BUS_STOP_N"))

p$diff <- "TRUE"
colnames(p)[colnames(p) == "AM_diff_unique"] <- "grid_id"
pp <- left_join(honeycomb_grid_sf, p)
ppp <- pp[!is.na(pp$diff), ]

mapA <- tm_shape(ppp) +
  tm_fill(fill = "Blues",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekday 6am-9am",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type = "8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha = 0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))

tmap_mode("view")
mapA

Here we see areas that underwent a change in distribution patterns. Most of the difference can be seen to be in the central business district where we expect to have higher crowd on weekdays as compared to the weekends.

Additional test

# hunan$Z.GDPPC <- scale(hunan$GDPPC) %>% 
#   as.vector 
# 
# 
# 
# nci2 <- moran.plot(hunan$Z.GDPPC, rswm_q,
#                    labels=as.character(hunan$County),
#                    xlab="z-GDPPC 2012", 
#                    ylab="Spatially Lag z-GDPPC 2012")

pacman::p_load(sf,tmap,spdep,tidyverse, dplyr, mapview, sfdep)